function [s_Yresdiff21, s_Ystatsdiff21, s_Xresdiff21] = ...
    TideDiff( s_Yres2, s_Yres1, s_Ystats2, s_Ystats1, s_Xresdiff21 );
%function [s_Yresdiff21, s_Ystatsdiff21, s_Xresdiff21] = ...
%    TideDiff( s_Yres2, s_Yres1, s_Ystats2, s_Ystats1, s_Xresdiff21 );
%
%Calculate the differences level2 - level1
%
%Output:
%s_Yresdiff21 = struct( ...
%        'signal', s_Yres2.signal(index) - s_Yres1.signal(index) ...
%    ,   'residual', s_Yres2.residual(index) - s_Yres1.residual(index) ...
%    ,   'tidepred', s_Yres2.tidepred(index) - s_Yres1.tidepred(index) ...
%    ,   'amplitude', s_Yres2.amplitude(index) - s_Yres1.amplitude(index) ...
%    ,   'phase', s_Yres2.phase(index) - s_Yres1.phase(index) ...
%);
%
%s_Ystatsdiff21 = struct( ...
%        'snr', s_Ystats2.snr ...
%    ,   'snr_crit', s_Ystats2.snr_crit ...
%    ,   'amperror', s_Ystats2.amperror + s_Ystats1.amperror ...
%    ,   'phaserror', s_Ystats2.phaserror + s_Ystats1.phaserror ...
%    ,   'rmsesignal', std(s_Yresdiff21.signal(index)) ...
%    ,   'rmseresidual', std(s_Yresdiff21.residual(index)) ...
%    ,   'rmsetidepred', std(s_Yresdiff21.tidepred(index)) ...
%    ,   'corrsignal', corrcoef(s_Ystats2.signal(index), s_Ystats1.signal(index)) ...
%    ,   'corrresidual', corrcoef(s_Ystats2.residual(index), s_Ystats1.residual(index)) ...
%    ,   'corrtidepred', corrcoef(s_Ystats2.tidepred(index), s_Ystats1.tidepred(index)) ...
%);

%Look for non-NaN indexes
index = find(~isnan(s_Yres2.signal));

s_Yresdiff21 = struct( ...
        'signal', s_Yres2.signal(index) - s_Yres1.signal(index) ...
    ,   'residual', s_Yres2.residual(index) - s_Yres1.residual(index) ...
    ,   'tidepred', s_Yres2.tidepred(index) - s_Yres1.tidepred(index) ...
    ,   'amplitude', s_Yres2.amplitude - s_Yres1.amplitude ...
    ,   'phase', s_Yres2.phase - s_Yres1.phase ...
);

s_Ystatsdiff21 = struct( ...
        'snr', s_Ystats2.snr ...
    ,   'snr_crit', s_Ystats2.snr_crit ...
    ,   'amperror', s_Ystats2.amperror + s_Ystats1.amperror ...
    ,   'phaserror', s_Ystats2.phaserror + s_Ystats1.phaserror ...
    ,   'rmsesignal', std(s_Yresdiff21.signal) ...
    ,   'rmseresidual', std(s_Yresdiff21.residual) ...
    ,   'rmsetidepred', std(s_Yresdiff21.tidepred) ...
    ,   'corrsignal', corrcoef(s_Yres2.signal(index), s_Yres1.signal(index)) ...
    ,   'corrresidual', corrcoef(s_Yres2.residual(index), s_Yres1.residual(index)) ...
    ,   'corrtidepred', corrcoef(s_Yres2.tidepred(index), s_Yres1.tidepred(index)) ...
);

s_Xresdiff21.time = s_Xresdiff21.time(index);

disp([  'The rmse and correlation of the original signals are ' ...
        num2str(s_Ystatsdiff21.rmsesignal) ...
        '(m) and ' ...
        num2str(s_Ystatsdiff21.corrsignal(1,2)) ...
     ]);
 
disp([  'The rmse and correlation of the signals tide-only component are ' ...
        num2str(s_Ystatsdiff21.rmsetidepred) ...
        '(m) and ' ...
        num2str(s_Ystatsdiff21.corrtidepred(1,2)) ...
     ]);
 
disp([  'The rmse and correlation of the tidefree signals are ' ...
        num2str(s_Ystatsdiff21.rmseresidual) ...
        '(m) and ' ...
        num2str(s_Ystatsdiff21.corrresidual(1,2)) ...
     ]);